Fitting waggle dance models to waggle dance data

This document contains code to fit the waggle dance model to each of our 20 different sites. The code here uses the wagglefit package, as well as some additional code stored in fit_data.R to help simplify things. Each site has its own section, with reused code. The final sections create map plots and summary statistic plots.

Optimising each site

# set waggle dance duration in seconds as foraging distance for analysis
data <- read.csv("data/FullHBForagingData.csv")

alldata <- data %>%
  filter(Year == 2017) %>% # remove pilot data conducted in 2016
  select(date, site, duration.seconds) %>%
  rename(foraging_distance = duration.seconds)

BEL

Provides a good fit on the data. All parameters look central and nicely covered.

target_site <- "BEL"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 1.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 0.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 8, -5.5, -1.8)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
BEL	collective	-181.4148	0.6475794	7.041850	0.8144976	0.0514862	0.1698127	5	372.8296	0.0980392	0.714
BEL	individual	-183.5011	1.0000000	4.495285	NA	0.0725148	NA	2	371.0021	0.0882353	0.798

all_sites[[target_site]]$fit

BFI

All parameters look central in the likelihood space and a nice fit is returned.

target_site <- "BFI"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 0.5)
bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 0.6)
bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.2, 3.5, -6.5, -2.)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
BFI	collective	-190.4990	0.0448114	0.000001	2.712405	0.1645321	0.2509062	5	390.9980	0.0873362	0.344
BFI	individual	-211.2967	1.0000000	9.470271	NA	0.1000000	NA	2	426.5934	0.1834061	0.002

all_sites[[target_site]]$fit

BLO

The individual fit isnt great and the Bs paramater is increasing up to the boundry, indicating it can only really take on a straight line / exponential fit. The parameters for the collective model are fairly central in the likelihood space, the fit looks very good.

target_site <- "BLO"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-10, 5)
br_bnds <- c(1.0e-10, 5)
as_bnds <- c(1.0e-5, 5)
ar_bnds <- c(1.0e-10, 1.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 5505)
as_bnds <- c(1.0e-12, 0.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.2, 9.5, -8., -2.2)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
BLO	collective	-221.0260	0.2085372	1.785832	0.570999	0.0677377	0.3591282	5	452.0519	0.0705882	0.768
BLO	individual	-237.3888	1.0000000	5505.000000	NA	0.0001222	NA	2	478.7776	0.1117647	0.215

all_sites[[target_site]]$fit

BUR

Bs and Br go in oposite directions. E.g. Bs approaches 0 whilst Br approaches an every higher number.

Fit looks good.

target_site <- "BUR"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 500)
br_bnds <- c(1.0e-6, 500)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 1.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 0.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 4, -6.8, -2.5)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords,
  xlims = c(0, 6)
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
BUR	collective	-118.3048	0.1116261	1e-06	499.0741	0.1513544	0.0275011	5	246.6096	0.0670732	0.829
BUR	individual	-143.7177	1.0000000	1e+01	NA	0.1000000	NA	2	291.4355	0.1768293	0.016

all_sites[[target_site]]$fit

CAD

The collective model roughly follows the individual model but is able to acheive a slightly improved fit to the shoulder, hence it provides a higher likelihood score. The AIC indicates this is overfitting, suggesting the individual model provides a more parsimonious explanation.

\(p\) does not approach 1 as one might expect, indicating there is overfitting and so a comparison with an individual model is required.

target_site <- "CAD"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 1.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 2)
as_bnds <- c(1.0e-12, 0.8)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 1.8, -6, -1.2)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
CAD	collective	-48.27822	0.164847	1e-06	0.425631	1.4569294	0.3860464	5	106.5564	0.0821918	0.949
CAD	individual	-49.95313	1.000000	1e-06	NA	0.3854731	NA	2	103.9063	0.1095890	0.745

all_sites[[target_site]]$fit

GIL

The collective model provides the most parsimonious and best fit.

target_site <- "GIL"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 1.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 0.8)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 2.1, -8., -2)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
GIL	collective	-74.54143	0.2614397	0.000001	1e-06	0.3342945	0.7547164	5	159.0829	0.1354167	0.053
GIL	individual	-102.30919	1.0000000	2.036505	NA	0.3450190	NA	2	208.6184	0.1979167	0.002

all_sites[[target_site]]$fit

HER

The collective model provides the best fit to the data, but the proportion of scouts is high.

target_site <- "HER"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 100)
br_bnds <- c(1.0e-6, 100)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 1.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 2)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 4, -6.5, -1.5)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
HER	collective	-136.9873	0.249628	1e-06	0.2155781	0.3956663	0.181939	5	283.9745	0.0674157	0.975
HER	individual	-138.9302	1.000000	1e-06	NA	0.1751111	NA	2	281.8604	0.1123596	0.565

all_sites[[target_site]]$fit

HHS

Collective model provides the best fit but falls under the tail and shoulder. The individual model strugles to find a good fit, probably due to the tail.

target_site <- "HHS"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 50)
br_bnds <- c(1.0e-6, 50)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 1.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 4, -6.5, -2.2)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
HHS	collective	-61.34799	0.0953466	0.000001	1e-06	0.1411228	0.5028644	5	132.6960	0.1333333	0.364
HHS	individual	-83.74915	1.0000000	7.095173	NA	0.1251512	NA	2	171.4983	0.2444444	0.006

all_sites[[target_site]]$fit

HOR

Collective model provides a very good fit, whilst the individual model fails to find much traction. The proportion of scouts goes very low (~3%) suggesting the majority of the colony are following a small number of scouting individuals.

target_site <- "HOR"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 1.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 3, -9, -2.5)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
HOR	collective	-40.50196	0.0557038	0.000001	0.5393193	0.2182387	0.6281992	5	91.00392	0.0533333	0.965
HOR	individual	-56.65192	1.0000000	7.774875	NA	0.2019358	NA	2	117.30385	0.1666667	0.032

all_sites[[target_site]]$fit

MAK

Again the collective model provides a good fit to the data, however, the individual model fits poorly, reducing to an exponential.

target_site <- "MAK"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 0.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 20)
as_bnds <- c(1.0e-12, 1.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# run individual model
individual_result <- fit_individual_model_to_data(data, bounds)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

# view individual model likelihood space to check bounds look ok
individual_result$llspace

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 4.2, -6.5, -2.5)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
MAK	collective	-84.51875	0.2093454	2.319365	0.0201563	0.1221729	0.4776196	5	179.0375	0.0909091	0.776
MAK	individual	-98.02457	1.0000000	9.168368	NA	0.0961027	NA	2	200.0491	0.1919192	0.044

all_sites[[target_site]]$fit

MEL

Collective model provides the best fit, however it misses a large section of the shoulder for the tail. The \(bs\) parameter wants to go to zero, however when let go bellow 1e-6 the behaviour becomes very erratic and the fit deteriorates.

target_site <- "MEL"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 0.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 3., -7, -2.1)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
MEL	collective	-115.3164	0.0806865	1e-06	1e-06	0.2128244	0.3377644	5	240.6327	0.1626016	0.088
MEL	individual	-135.7308	1.0000000	1e-06	NA	0.2520904	NA	2	275.4617	0.2439024	0.001

all_sites[[target_site]]$fit

MPA

Collective provides the best fit.

target_site <- "MPA"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 0.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 4., -7.2, -2.1)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
MPA	collective	-181.8729	0.5669623	0.5783197	1e-06	0.1853651	0.4199909	5	373.7458	0.0666667	0.903
MPA	individual	-188.0378	1.0000000	1.4738536	NA	0.1866656	NA	2	380.0756	0.1000000	0.397

all_sites[[target_site]]$fit

ROT

Colletive provides the best fit.

target_site <- "ROT"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 0.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 2., -5.5, -2.1)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
ROT	collective	-138.1293	0.309561	1e-06	1e-06	0.3151219	0.4832754	5	286.2585	0.1546392	0.018
ROT	individual	-153.3550	1.000000	1e-06	NA	0.3538107	NA	2	310.7100	0.2010309	0.002

all_sites[[target_site]]$fit

SAU

Collective provides the best fit.

target_site <- "SAU"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 0.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 2.5, -6.5, -2.1)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
SAU	collective	-108.1850	0.3096488	0.000001	10	0.2620339	0.1638241	5	226.3700	0.0606061	0.946
SAU	individual	-115.5025	1.0000000	1.488528	NA	0.2710389	NA	2	235.0051	0.1287879	0.197

all_sites[[target_site]]$fit

SOM

Collective provides the best fit but misses a large section of the shoulder.

target_site <- "SOM"

# subset data for target site
data <- alldata %>%
  filter(site == target_site) %>%
  filter(foraging_distance < 6) # remove outlier distance

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 0.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 2., -5.5, -2.1)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
SOM	collective	-70.53597	0.0020792	1e+01	10	1.5000000	0.1486462	5	151.0719	0.1071429	0.492
SOM	individual	-76.67133	1.0000000	1e-06	NA	0.3772311	NA	2	157.3427	0.1607143	0.087

all_sites[[target_site]]$fit

SRA

Collective provides the best fit.

target_site <- "SRA"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 0.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 3., -7.5, -2.1)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
SRA	collective	-123.3686	0.2424498	0.000001	0.205498	0.2128967	0.492155	5	256.7373	0.0519481	0.982
SRA	individual	-138.1142	1.0000000	3.002015	NA	0.2133652	NA	2	280.2284	0.1428571	0.082

all_sites[[target_site]]$fit

STU

Individual provides the best fit. Again, the proportion of scouts does not move towards 1. I.e. the collective model fails to reduce to the individual in the MLE fit.

target_site <- "STU"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 0.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 4, -6.5, -1.5)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
STU	collective	-155.7499	0.7988692	0.4228714	0.3278325	0.1638306	0.5	5	321.4997	0.0454545	0.999
STU	individual	-156.5403	1.0000000	1.3465643	NA	0.1612383	NA	2	317.0806	0.0545455	0.997

all_sites[[target_site]]$fit

SWP

Collective provides the best fit.

target_site <- "SWP"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 15)
br_bnds <- c(1.0e-6, 15)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 0.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 2., -6, -1.5)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
SWP	collective	-40.42382	0	2.706000	1.274713	0.2695819	0.3841755	5	90.84765	0.1000000	0.731
SWP	individual	-46.44077	1	0.000001	NA	0.4236173	NA	2	96.88154	0.1444444	0.260

all_sites[[target_site]]$fit

YAL

Collective provides the best fit.

target_site <- "YAL"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 0.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.1)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 4., -8, -2.)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
YAL	collective	-205.7605	0.0786728	0.0000010	0.3007459	0.1812311	0.2730218	5	421.5211	0.0439560	0.991
YAL	individual	-222.4328	1.0000000	0.2168113	NA	0.2128333	NA	2	448.8656	0.0989011	0.318

all_sites[[target_site]]$fit

ZSL

Collective provides the best fit.

target_site <- "ZSL"

# subset data for target site
data <- alldata %>%
  filter(site == target_site)

# set up bounds for the collective model
p_bnds <- c(0, 1.0)
bs_bnds <- c(1.0e-6, 10)
br_bnds <- c(1.0e-6, 10)
as_bnds <- c(1.0e-12, 1.5)
ar_bnds <- c(1.0e-12, 0.5)
collective_bounds <- rbind(
  p_bnds, bs_bnds,
  br_bnds, as_bnds,
  ar_bnds
)

# set up bounds for the individual model
bs_bnds <- c(1.0e-6, 50)
as_bnds <- c(1.0e-12, 1.3)
individual_bounds <- rbind(
  bs_bnds, as_bnds
)

# set coordinates for histogram subplot
subplot_coords <- c(0.1, 1.8, -8, -2)

all_sites[[target_site]] <- run_wagglefit_analysis(
  target_site, data, collective_bounds, individual_bounds, subplot_coords
)

#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"
#> [1] "Itteration 1"
#> [1] "Itteration 2"
#> [1] "Itteration 3"
#> [1] "Itteration 4"
#> [1] "Itteration 5"
#> [1] "Itteration 6"
#> [1] "Itteration 7"
#> [1] "Itteration 8"
#> [1] "Itteration 9"
#> [1] "Itteration 10"

all_sites[[target_site]]$fit_result %>%
  kbl() %>%
  kable_classic(full_width = F)

site	model	loglikelihood	p	bs	br	as	ar	k	AIC	ks_statistic	ks_pvalue
ZSL	collective	-109.9979	0.2571773	0.0000010	0.8022074	0.3868084	0.5	5	229.9958	0.0361991	1.000
ZSL	individual	-119.9481	1.0000000	0.3918133	NA	0.4194463	NA	2	243.8962	0.0950226	0.275

all_sites[[target_site]]$fit

Overall findings

# save the analysis results for all sites
# save(all_sites, file = "results/site_fit_results.Rdata")

# group all site results together
df <- map(all_sites, 1) %>%
  bind_rows()

# save results
# saveRDS(df, file = "results/site_fit_results.Rda")

# all_sites$BFI$fit <- all_sites$BFI$fit +
#   ylim(0, 6) +
#   scale_x_continuous(breaks = seq(0, 6, by = 2))

# AIC plot
aic_plot <- df %>%
  group_by(site) %>%
  slice(which.min(AIC)) %>%
  select(model) %>%
  group_by(model) %>%
  summarise(lowest_AIC = n()) %>%
  ggplot(aes(x = model, y = lowest_AIC)) +
  geom_bar(stat = "identity") +
  labs(x = "Model", y = "Count") +
  scale_y_continuous(breaks = seq(0, 20, by = 2)) +
  theme(
    text = element_text(size = 42)
  )

ggsave(
  plot = aic_plot,
  filename = "results/figures/AIC_plot.png",
  width = 90,
  height = 110,
  units = "mm",
  dpi = 300
)

# ks plot
ks_plot_dist <- df %>%
  ggplot(aes(x = ks_pvalue)) +
  geom_histogram(bins = 10, binwidth = 0.1, col = "white") +
  geom_vline(xintercept = 0.05, color = "red", linetype = "dashed") +
  scale_y_continuous(breaks = seq(0, 12, by = 2)) +
  labs(x = "KS P value") +
  facet_wrap(~model, nrow = 3) +
  theme(
    text = element_text(size = 42),
    strip.background = element_blank()
  )


ggsave(
  plot = ks_plot_dist,
  filename = "results/figures/sites_ks.png",
  width = 86,
  height = 180,
  units = "mm",
  dpi = 300
)

ggsave(
  plot = all_sites$STU$fit,
  filename = "results/figures/STU.png",
  width = 90,
  height = 110,
  units = "mm",
  dpi = 300
)

ggsave(
  plot = all_sites$BFI$fit,
  filename = "results/figures/BUR.png",
  width = 90,
  height = 110,
  units = "mm",
  dpi = 300
)

Map plots

Code to make the individual mal plots. These figures are created standalone but are combined in to facets manually in an image processor.

library(ggplot2)
library(ggrepel)
library(gridExtra)
library(ggsn)
library(sf)
library(rworldmap)
library(ggspatial)
library(rnaturalearth)
library(rnaturalearthdata)


full_data_path <- "data/FullHBForagingData.csv"
data_raw <- tibble(read.csv(full_data_path))

map_data <- data_raw %>%
  select(
    site, lat, lon
  ) %>%
  distinct()

map_data <- df %>%
  group_by(site) %>%
  slice(which.min(AIC)) %>%
  select(site, model) %>%
  left_join(map_data, on = "site") %>%
  mutate(col = ifelse(site %in% c("STU", "BUR"), "1", "0"))

locations <- st_as_sf(
  map_data,
  coords = c("lon", "lat"), crs = 4326
)

# Extract selected sites for figure
selected_sites <- filter(map_data, site %in% c("STU", "BUR")) %>%
  mutate(
    label = ifelse(site == "STU", "c", "d")
  )

points_area <- st_bbox(locations)

worldmap <- ne_countries(scale = "large", returnclass = "sf")

# london area
london <- st_read(
  "shapefiles/London_Ward.shp"
)

#> Reading layer `London_Ward' from data source 
#>   `/home/joe/Documents/phd/HoneybeeResearch/wagglefit/analysis/shapefiles/London_Ward.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 649 features and 0 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 503568.2 ymin: 155850.8 xmax: 561957.5 ymax: 200933.9
#> Projected CRS: OSGB 1936 / British National Grid

london <- st_transform(
  london,
  CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
)

inset <- ggplot() +
  geom_sf(
    data = worldmap,
    fill = "grey90",
    color = "#4b4949d0"
  ) +
  geom_sf(
    data = london,
    fill = "#ced0cffc",
    lwd = 0
  ) +
  coord_sf(
    xlim = c(points_area[[1]] - 0.1, points_area[[3]] + 0.1),
    ylim = c(points_area[[2]] - 0.1, points_area[[4]] + 0.1)
  ) +
  geom_point(
    data = map_data,
    aes(x = lon, y = lat, shape = model, colour = model), size = 1.5
  ) +
  geom_text(
    data = selected_sites, aes(x = lon, y = lat, label = label),
    nudge_x = c(0, -0.07), nudge_y = c(0.05, 0.05), size = 8
  ) +
  geom_point(
    data = selected_sites, aes(x = lon, y = lat), size = 1.5, shape = 4
  ) +
  annotation_north_arrow(
    location = "tr", which_north = "true",
    style = north_arrow_fancy_orienteering,
    height = unit(15, "mm"),
    width = unit(15, "mm"),
    text_cex = 1.5
  ) +
  annotation_scale(
    location = "br",
    text_cex = 1.5
  ) +
  scale_shape_manual(values = c(1, 2)) +
  scale_colour_manual(values = c("black", "red")) +
  theme_nothing() +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(x = NULL, y = NULL) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = .5)
  )

inset

# make full plot
base <- ggplot(data = worldmap) +
  geom_sf(
    fill = "#c4cfc8",
    color = "#4b4949d0",
    lwd = 0.2
  ) +
  coord_sf(
    xlim = c(-11, 3),
    ylim = c(49.5, 60)
  ) +
  # geom_point(data = map_data, aes(x = lon, y = lat), size = 0.2) +
  geom_rect(
    aes(
      xmin = points_area$xmin[[1]] - 0.1, xmax = points_area$xmax[[1]] + 0.1,
      ymin = points_area$ymin[[1]] - 0.1, ymax = points_area$ymax[[1]] + 0.1
    ),
    fill = NA,
    colour = "black",
    size = .02
  ) +
  theme_nothing() +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(x = NULL, y = NULL) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = .5)
  )

base

merge_plot <- base +
  annotation_custom(
    ggplotGrob(inset),
    xmin = 1,
    xmax = 13,
    ymin = 52.5,
    ymax = 60
  )

merge_plot

ggsave(
  plot = merge_plot,
  filename = "results/figures/site_map.png",
  width = 90,
  height = 110,
  units = "mm",
  dpi = 300
)


merge_plot_2 <- inset +
  annotation_custom(
    ggplotGrob(base),
    xmin = -3.45,
    xmax = 2.2,
    ymin = 51.01,
    ymax = 51.45
  )

merge_plot_2

sites_model_plot <- plot_grid(
  merge_plot_2, ks_plot_dist, all_sites$STU$fit, all_sites$BUR$fit,
  labels = c("A", "B", "C", "D"), label_size = 24
)

ggsave(
  plot = sites_model_plot,
  filename = "results/figures/sites_model_plot.png",
  width = 183,
  height = 190,
  units = "mm",
  dpi = 300
)

ggsave(
  plot = sites_model_plot,
  filename = "results/figures/sites_model_plot.svg",
  width = 183,
  height = 190,
  units = "mm",
  dpi = 300
)

ggsave(
  plot = sites_model_plot,
  filename = "results/figures/sites_model_plot.pdf",
  width = 183,
  height = 190,
  units = "mm",
  dpi = 300
)

ggsave(
  plot = sites_model_plot,
  filename = "results/figures/sites_model_plot.bmp",
  width = 183,
  height = 190,
  units = "mm",
  dpi = 600
)

model_sites_fits <- plot_grid(
  merge_plot_2, ks_plot_dist,
  labels = c("A", "B"), label_size = 24
)

ggsave(
  plot = model_sites_fits,
  filename = "results/figures/model_sites_fits.png",
  width = 183,
  height = 190,
  units = "mm",
  dpi = 300
)

stu_zsl_fit <- plot_grid(
  all_sites$STU$fit, all_sites$BUR$fit,
  labels = c("A", "B"), label_size = 24
)

ggsave(
  plot = stu_zsl_fit,
  filename = "results/figures/stu_zsl_fit.svg",
  width = 183,
  height = 190,
  units = "mm",
  dpi = 100
)

# other map plots
inset <- ggplot() +
  geom_sf(
    data = worldmap,
    fill = "grey90",
    color = "#4b4949d0"
  ) +
  geom_sf(
    data = london,
    fill = "#ced0cffc",
    lwd = 0
  ) +
  coord_sf(
    xlim = c(points_area[[1]] - 0.1, points_area[[3]] + 0.1),
    ylim = c(points_area[[2]] - 0.1, points_area[[4]] + 0.1)
  ) +
  geom_point(
    data = map_data,
    aes(x = lon, y = lat, shape = model, colour = model), size = 1.5
  ) +
  annotation_north_arrow(
    location = "tr", which_north = "true",
    style = north_arrow_fancy_orienteering,
    height = unit(10, "mm"),
    width = unit(10, "mm")
  ) +
  annotation_scale(
    location = "br",
    text_cex = 3
  ) +
  scale_shape_manual(values = c(1, 2)) +
  scale_colour_manual(values = c("black", "red")) +
  theme_nothing() +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(x = NULL, y = NULL) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = .5),
    text = element_text(size = 42)
  )

inset

# make full plot
base <- ggplot(data = worldmap) +
  geom_sf(
    fill = "#c4cfc8",
    color = "#4b4949d0",
    lwd = 0.2
  ) +
  coord_sf(
    xlim = c(-11, 15),
    ylim = c(49.5, 60)
  ) +
  geom_point(data = map_data, aes(x = lon, y = lat), size = 0.2) +
  geom_rect(
    aes(
      xmin = points_area$xmin[[1]] - 0.1, xmax = points_area$xmax[[1]] + 0.1,
      ymin = points_area$ymin[[1]] - 0.1, ymax = points_area$ymax[[1]] + 0.1
    ),
    fill = NA,
    colour = "black",
    size = .02
  ) +
  theme_nothing() +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(x = NULL, y = NULL) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = .5)
  )

merge_plot <- base +
  annotation_custom(
    ggplotGrob(inset),
    xmin = 1,
    xmax = 13,
    ymin = 52.5,
    ymax = 60
  )

merge_plot

ggsave(
  plot = merge_plot,
  filename = "results/figures/site_map.png",
  width = 183,
  height = 190,
  units = "mm",
  dpi = 300
)


sites_model_plot <- plot_grid(
  all_sites$STU$fit, all_sites$BUR$fit, aic_plot, ks_plot_dist,
  labels = c("A", "B", "C", "D"), label_size = 42
)

sites_model_plot

ggsave(
  plot = sites_model_plot,
  filename = "results/figures/results_model_plot.png",
  width = 183,
  height = 190,
  units = "mm",
  dpi = 300
)